{
 "metadata": {
  "name": ""
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "heading",
     "level": 1,
     "metadata": {},
     "source": [
      "Teste da nova clase Gas"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Se ha escrito um novo programa para a simula\u00e7\u00e3o do fluxo de calor num gas bidimensional de Lennard-Jones. Porque quero provar o seguintes casos:\n",
      "- Lennard-Jones SF\n",
      "- Armonico\n"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "$$ \\Phi^{LJ}(r) = -4 \\varepsilon \\left( \\left(\\frac{\\sigma}{r}\\right)^{12} - \\left(\\frac{\\sigma}{r}\\right)^6 \\right) $$"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "$$\\frac{\\partial }{\\partial r} \\Phi^{LJ}(r) = -24 \\varepsilon \\left(2 \\left(\\frac{\\sigma}{r}\\right)^{12} - \\left(\\frac{\\sigma}{r}\\right)^6 \\right)\\frac{1}{r} $$"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "$$F^{LJ}_x(r) = -\\frac{\\partial }{\\partial r} \\Phi^{LJ}(r) \\frac{\\partial r}{\\partial x}$$"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "$$ \\Phi^{SF}(r) = \\Phi^{LJ}(r) - \\Phi^{LJ}(rc) - \\frac{\\partial}{\\partial r} \\Phi^{LJ}(r)\\Big|_{r=rc} (\\frac{r^2 - rc^2}{2rc}) $$"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "$$F^{SF}_x(r) = -\\frac{\\partial }{\\partial r} \\Phi^{SF}(r) \\frac{\\partial r}{\\partial x}$$"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "$$ r = \\sqrt{x^2 + y^2 +z^2} $$"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "$$ \\frac{\\partial r}{\\partial x} = \\frac{1}{r} x $$"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "$$\n",
      "F^{SF}_x(r) = \n",
      "\\left(\n",
      "-\n",
      "\\frac{\\partial }{\\partial r} \\Phi^{LJ}(r) \n",
      "+ \n",
      "\\frac{\\partial }{\\partial r} \\Phi^{LJ}(r)\\Big|_{r=rc} \\left(\\frac{r}{rc} \\right)\n",
      "\\right)\n",
      "\\frac{x}{r}\n",
      "$$"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "$$\n",
      "F^{SF}_x(r) = \n",
      "\\left(\n",
      " -\\left(\\frac{1}{r} \\right)\n",
      "\\frac{\\partial }{\\partial r} \\Phi^{LJ}(r)\n",
      "+ \n",
      " \\left(\\frac{1}{rc} \\right) \\frac{\\partial }{\\partial r} \\Phi^{LJ}(r)\\Big|_{r=rc}\n",
      "\\right) x\n",
      "$$"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "$$\n",
      "DrPhiLJOverr =\n",
      "\\left(\\frac{1}{r} \\right)\n",
      "\\frac{\\partial }{\\partial r} \\Phi^{LJ}(r) \n",
      "$$"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "$$\n",
      "F^{LJ}_x(r) = \n",
      "\\left(\n",
      "24 \\varepsilon \\left(2 \\left(\\frac{\\sigma}{r}\\right)^{12} - \\left(\\frac{\\sigma}{r}\\right)^6 \\right)\\frac{1}{r}\n",
      "+ \n",
      "\\frac{\\partial }{\\partial r} \\Phi^{LJ}(r)\\Big|_{r=rc} \\left(\\frac{r}{rc} \\right)\n",
      "\\right)\n",
      "\\frac{x}{r}\n",
      "$$"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "$$\n",
      "F^{LJ}_x(r) = \n",
      "\\left(\n",
      "24 \\varepsilon \\left(2 \\left(\\frac{\\sigma}{r}\\right)^{12} - \\left(\\frac{\\sigma}{r}\\right)^6 \\right) \\frac{x}{r^2}\n",
      "+ \n",
      "\\frac{\\partial }{\\partial r} \\Phi^{LJ}(r)\\Big|_{r=rc} \\left(\\frac{x}{rc} \\right)\n",
      "\\right)\n",
      "$$"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%run Gas.py"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "gas_old = Gas(9)\n",
      "gas_new = Gas(9)\n",
      "\n",
      "gas_ini = copy.deepcopy(gas_old)\n",
      "\n",
      "gamma  = 1.0\n",
      "Tleft  = 1.0\n",
      "Tright = 1.5\n",
      "Wrange = 1.0\n",
      "dt = 0.001\n",
      "\n",
      "total_time = 100\n",
      "\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "gas_ini.vec_v"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 3,
       "text": [
        "array([[-1.75989115,  0.        ,  0.        ],\n",
        "       [-2.49976415,  0.        ,  0.        ],\n",
        "       [ 0.97637169,  0.        ,  0.        ],\n",
        "       [ 1.1664947 ,  0.        ,  0.        ],\n",
        "       [ 2.47050558,  0.        ,  0.        ],\n",
        "       [ 0.98957054,  0.        ,  0.        ],\n",
        "       [ 0.5135044 ,  0.        ,  0.        ],\n",
        "       [-0.61740711,  0.        ,  0.        ],\n",
        "       [-0.65602709,  0.        ,  0.        ]])"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "total_time =10000\n",
      "enerDensTemp = np.zeros(total_time*gas_old.N).reshape(total_time,gas_old.N)\n",
      "fluxDensTemp = np.zeros(total_time*gas_old.N).reshape(total_time,gas_old.N)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 26
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Simulation para um tempo total_time"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "######\n",
      "for i in range(total_time):\n",
      "    gas_old.calculateDeltas()\n",
      "    gas_old.calculateForces()\n",
      "    gas_old.leftHB(gamma, Tleft, Wrange)\n",
      "    gas_old.rightHB(gamma, Tright, Wrange)\n",
      "    \n",
      "    gas_new.stepPosition(gas_old,dt)\n",
      "    gas_new.inBox()\n",
      "    \n",
      "    gas_new.calculateDeltas()\n",
      "    gas_new.calculateForces()\n",
      "    gas_new.stepVelocity(gas_old, dt)\n",
      "    \n",
      "    gas_new.calculateHalfV()\n",
      "    gas_new.calculateEnerDensity()  # h_n\n",
      "    gas_new.calculateFluxDensity()  # j_n\n",
      "    \n",
      "    #gas_new.calculateEpot()\n",
      "    #gas_new.calculateEkin()\n",
      "    #gas_new.calculateFlux()\n",
      "    \n",
      "    #save the values for the next step\n",
      "    gas_old = copy.deepcopy(gas_new)\n",
      "    gas_old.calculateEkin_n()\n",
      "    enerDensTemp[i] = gas_old.Ekin_n\n",
      "    gas_old.calculateFluxDensity()\n",
      "    fluxDensTemp[i] = gas_old.fluxDensity[:,0]\n",
      "    #plt.plot(range(len(gas_old.Ekin_n)), gas_old.Ekin_n)\n",
      "\n",
      "######\n",
      "gas_end = copy.deepcopy(gas_new)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 28
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "fig, ax = plt.subplots()\n",
      "plt.title(\"$ j_{x_n} $\")\n",
      "plt.ylabel(\"$t$ steps\")\n",
      "plt.xlabel(\"$n$\")\n",
      "#plt.colorbar()\n",
      "#heatmap = ax.pcolor(enerDensTemp[:,1:gas_old.N-2], cmap='RdBu')#plt.cm.Blues)\n",
      "heatmap = ax.pcolor(fluxDensTemp[:,1:gas_old.N-2], cmap='RdBu')#plt.cm.Blues)\n",
      "fig.colorbar(heatmap)\n",
      "#plt.imshow(enerDensTemp[:,1:gas_old.N-2])\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 32
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "gas_new.vec_v"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 35,
       "text": [
        "array([[-0.7034278 ,  0.        ,  0.        ],\n",
        "       [ 0.0400347 ,  0.        ,  0.        ],\n",
        "       [ 0.05377395,  0.        ,  0.        ],\n",
        "       [ 0.05917797,  0.        ,  0.        ],\n",
        "       [ 0.11446098,  0.        ,  0.        ],\n",
        "       [ 0.14731351,  0.        ,  0.        ],\n",
        "       [ 0.10447626,  0.        ,  0.        ],\n",
        "       [ 0.07826526,  0.        ,  0.        ],\n",
        "       [ 0.00104247,  0.        ,  0.        ]])"
       ]
      }
     ],
     "prompt_number": 35
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#plt.plot(range(len(gas_old.enerDensity)), gas_old.enerDensity)\n",
      "plt.plot(range(len(gas_old.Ekin_n)), gas_old.Ekin_n)\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 36
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "gas_ini.plot()\n",
      "gas_end.plot()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 11
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "fluxDensTemp[-1000]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 45,
       "text": [
        "array([-0.00351499,  0.01019262,  0.00835847,  0.0078035 ,  0.00535785,\n",
        "        0.00140925, -0.00361943, -0.02243463, -0.02738639])"
       ]
      }
     ],
     "prompt_number": 45
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "gas_end.fluxDensity"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 42,
       "text": [
        "array([[-0.02359147,  0.        ,  0.        ],\n",
        "       [-0.00103691, -0.        , -0.        ],\n",
        "       [-0.00209368, -0.        , -0.        ],\n",
        "       [-0.0030872 , -0.        , -0.        ],\n",
        "       [-0.0025165 , -0.        , -0.        ],\n",
        "       [-0.00565505, -0.        , -0.        ],\n",
        "       [-0.01016107,  0.        ,  0.        ],\n",
        "       [-0.01027885,  0.        ,  0.        ],\n",
        "       [-0.00300785, -0.        , -0.        ]])"
       ]
      }
     ],
     "prompt_number": 42
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "gas_ini.vec_r"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 9,
       "text": [
        "array([[-4. , -0.5,  0. ],\n",
        "       [-3. , -0.5,  0. ],\n",
        "       [-2. , -0.5,  0. ],\n",
        "       [-1. , -0.5,  0. ],\n",
        "       [ 0. , -0.5,  0. ],\n",
        "       [ 1. , -0.5,  0. ],\n",
        "       [ 2. , -0.5,  0. ],\n",
        "       [ 3. , -0.5,  0. ]])"
       ]
      }
     ],
     "prompt_number": 9
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "gas_ini.stepPosition(gas_ini,dt)\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 10
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print gas_ini.vec_v[:,0], gas_end.vec_v[:,0]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[-1.75989115 -2.49976415  0.97637169  1.1664947   2.47050558  0.98957054\n",
        "  0.5135044  -0.61740711 -0.65602709] [-0.7034278   0.0400347   0.05377395  0.05917797  0.11446098  0.14731351\n",
        "  0.10447626  0.07826526  0.00104247]\n"
       ]
      }
     ],
     "prompt_number": 52
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "gas_ini.calculateEkin_n()\n",
      "gas_end.calculateEkin_n()\n",
      "print gas_ini.Ekin_n.sum(), gas_end.Ekin_n.sum()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "9.90897338169 0.277325762361\n"
       ]
      }
     ],
     "prompt_number": 56
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "enerDensTemp.sum(axis=1)\n",
      "plt.plot(np.array(range(total_time))*dt,enerDensTemp.sum(axis=1))\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 71
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "gas_ini.force"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 8,
       "text": [
        "array([[ 0.,  0.,  0.],\n",
        "       [ 0.,  0.,  0.],\n",
        "       [ 0.,  0.,  0.],\n",
        "       [ 0.,  0.,  0.],\n",
        "       [ 0.,  0.,  0.],\n",
        "       [ 0.,  0.,  0.],\n",
        "       [ 0.,  0.,  0.],\n",
        "       [ 0.,  0.,  0.]])"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}